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ABSTRACT 

The  United  States  relies  heavily  on  its  space  infrastructure  for  a  vast  number  of  applications,  including 
communication,  navigation,  banking,  national  security,  and  research.  However,  NASA  predicts  that  between  now 
and  2030  orbital  collisions  will  become  increasingly  frequent  and  could  reach  a  runaway  environment.  This 
devastating  scenario,  also  known  as  the  Kessler  Syndrome,  has  the  potential  to  eventually  destroy  our  assets  in  near- 
Earth  space  and  result  in  a  debris  cloud  that  could  make  space  itself  inaccessible.  Preventing  the  Kessler  Syndrome 
requires,  in  addition  to  an  object  removal  technique,  a  groundbreaking  new  orbital  dynamics  framework  that 
combines  a  comprehensive  physics-based  model  of  atmospheric  drag  with  an  accurate  uncertainty  quantification  of 
orbital  predictions.  The  IMPACT  project  (Integrated  Modeling  of  Perturbations  in  Atmospheres  for  Conjunction 
Tracking),  funded  with  over  $5  Million  by  the  Los  Alamos  Laboratory  Directed  Research  and  Development  office, 
has  the  goal  to  develop  such  an  integrated  system  of  atmospheric  drag  modeling,  orbit  propagation,  and  conjunction 
analysis  with  detailed  uncertainty  quantification  to  address  the  space  debris  and  collision  avoidance  problem.  We 
discuss  the  components  and  capabilities  of  the  IMPACT  framework  and  show  a  short  demonstration  of  modeling 
interface  and  resulting  3D  visualizations. 


Project  PI:  Josef  Koller,  jkoller@lan.gov 


1.  INTRODUCTION 


The  IMPACT  (Figure  1)  project,  funded  by  an  internal  research  and  development  process  at  Los  Alamos  National 
Laboratory  (LANL),  has  the  goal  to  develop  an  integrated  modeling  system  for  addressing  current  needs  in  space 
debris  and  conjunction  analysis  for  resident  space  objects  (RSO)  in  low-Earth  orbit.  Now  with  almost  three  years 
into  the  project,  we  have  developed  an  integrated  solution  combining  physics-based  density  modeling  of  the  upper 
atmosphere  between  120-700  km  altitude,  satellite  drag  forecasting  for  quiet  and  disturbed  geomagnetic  conditions, 
and  conjunction  analysis  with  non-Gaussian  uncertainty  quantification.  We  are  employing  several  novel  approaches 
including 

•  a  Raven-class  observational  facility  for  ground-based  observations  of  space  objects 

•  physics-based  density  modeling  and  data  assimilation  for  integrating  density  measurements  with  model 
forecasts  similar  to  terrestrial  weather  prediction 

•  drag  coefficient  modeling  with  Test  Particle  and  Direct  Simulation  Monte  Carlo  methods  to  accurately 
account  for  changes  in  density,  chemical  composition,  temperature  and  their  effect  on  the  drag  coefficient 
and  orbit  propagation  during  geomagnetic  storms 

•  atmospheric  density  reconstruction  using  Satellite  Orbit  Tomography  based  on  X-ray  Computed 
Tomography  from  the  medical  imaging  field 

•  collision  statistics  and  uncertainty  quantification  using  full  Monte  Carlo  sampling  and  importance  sampling 
for  collisional  probabilities  and  allowing  for  non-Gaussian  probability  distributions 

•  machine  learning  approach  to  enable  a  coupling  between  solar  drivers  and  the  upper  atmosphere  resulting 
in  tremendous  computational  efficiency  while  preserving  modeling  accuracy 

m  m  i 

Integrated  Modeling  of  Perturbations  in  Atmospheres  for  Conjunction  Tracking 

Figure  1  IMPACT  logo 

In  the  following  subsection,  we  will  describe  in  details  some  of  the  technological  developments  as  part  of  the 
IMPACT  project. 


2.  GROUND  BASED  OBSERVATIONS 

LANL  is  using  a  Raven-class  telescope  (0.35  m  aperture  C14  on  a  Paramount  ME  mount)  to  track  satellites.  This 
observational  facility  is  located  at  2650  m  altitude  under  dark  skies  about  an  hour  from  Los  Alamos,  NM.  We 
typically  take  several  nights  of  observation  around  each  New  Moon,  weather  and  wildfire  permitting.  This  provides 
both  good  metric  accuracy  for  orbit  determination  and  lightcurves  for  object  characterization.  The  objects  we 
observe  for  IMPACT  are  primarily: 

a)  Satellites  with  accurately  known  orbits,  which  we  use  to  characterize  the  accuracy  of  our  metric  determinations. 
These  satellites  include  the  GPS  and  WAAS  constellations,  geodesy  satellites,  and  satellites  carrying  on-board  GPS 
receivers. 

b)  Cubesats,  many  of  which  are  similar  in  tenns  of  mass,  volume,  shape,  and  surface  materials  and  hence  are  good 
test  particles  for  looking  at  atmosphere-induced  variation  in  drag.  We  also  observe  other  LEO  satellites  that  have 
size/shape/mass  information  available. 

c)  Objects  in  highly  eccentric  orbits  that  dip  into  the  atmosphere  at  perigee.  The  drag  perturbation  over  an  orbit  for 
one  of  these  objects  is  dominated  by  the  atmosphere  at  the  geographic  location  and  altitude  of  the  perigee  pass. 
Rocket  bodies  in  GTO  are  numerous  and  particularly  useful,  since  they  tend  have  known  physical  characteristics. 
We  produce  light-curves  from  our  observations  in  order  to  derive  the  spin-states  of  the  objects,  which  affect  the 
average  frontal  area  around  the  time  of  perigee. 


We  have  also  observed  debris  from  satellite  collisions  and  the  Breeze-M  explosion.  These  observations  include  over 
300  different  RSOs,  including  26  GPS  satellites  and  76  rocket  bodies. 

3.  PHYSICS-BASED  DENSITY  FORECAST  MODELING 


A  key  element  in  improving  satellite  orbital  predictions  is  the  correct  specification  of  the  ionosphere-thermosphere 
environment,  given  that  atmospheric  density  exerts  significant  drag  over  satellites.  There  are  a  number  of  models 
that  can  estimate  the  composition  and  density  of  the  ionosphere-thermosphere,  from  empirical  models  to  physics 
based  models.  Empirical  models,  such  as  the  Mass  Spectrometer  and  Incoherent  Scatter  (MSIS)  [1]  model,  can 
provide  an  accurate  estimation  of  current  or  past  ionosphere-thermosphere  density,  based  on  a  number  of 
observations.  Unfortunately,  these  types  of  models  do  not  have  predictive  capabilities,  that  is,  they  provide  a  good 
nowcast  but  are  often  not  suitable  for  a  forecast.  The  IMPACT  project  uses  a  physics-based  model,  which  has  the 
potential  to  estimate  a  forecast  of  the  ionosphere-thermosphere  since  they  include  the  relevant  physical  behavior  of 
the  system.  In  particular,  the  project  uses  the  Global  Ionosphere-Thermosphere  Model  (GITM)  [2]  but  the  IMPACT 
framework  is  in  able  to  switch  between  different  models  for  comparative  studies. 

3.1  Ionosphere-Thermosphere  Models  and  Assimilation  Method 

GITM  is  a  physics-based  three-dimensional  model  that  solves  the  full  Navier-Stokes  equations  for  density, 
velocity,  and  temperature  for  a  number  of  neutral  and  charged  components.  To  account  for  solar  activity,  GITM  uses 
at  the  moment  the  FI 0.7  solar  flux,  hemispheric  power  index  (HPI)  (which  is  derived  from  the  3 -hour  Kp), 
interplanetary  magnetic  field  (IMF)  data  and  solar  wind  velocity.  GITM  inherently  allows  for  non-hydrostatic 
solutions  to  develop  which  allows  for  realistic  dynamics  in  the  auroral  zones  [2].  As  with  many  of  the  physics-based 
models,  GITM  includes  a  number  of  assumptions  and  physical  representation  of  the  ionosphere-thermosphere  that 
might  not  be  accurate  and  introduce  errors  into  the  estimation  of  the  density.  This  severely  affects  the  quality  of  a 
density  forecast. 

In  order  to  improve  the  predictability  of  GITM,  and  provide  a  good  forecast  of  the  ionosphere- 
thermosphere,  we  implement  a  data  assimilation  system  based  on  the  ensemble  Kalman  filter  (EnKF)  [3,4].  The 
EnKF  uses  an  ensemble  of  model  simulations  to  approximate  the  probability  distribution  of  the  model,  as  well  as  the 
covariance  matrix.  The  main  advantages  of  the  EnKF  are  the  ease  of  implementation  and  the  computational 
efficiency  for  non-linear  models.  In  particular  we  use  the  localized  ensemble  transform  Kalman  filter  (LETKF)  [5], 
which  is  a  localized  version  of  the  EnKF.  The  LETKF  assimilates  by  local  volume  centered  at  each  grid-point 
variable,  where  the  area  of  the  local  volume  depends  on  model  dynamics  and  assumptions  of  correlations  between 
model  variables.  Given  the  local  nature  of  the  LETKF,  the  algorithm  is  highly  parallel  since  all  grid-point  variables 
can  be  assimilated  simultaneously. 


3.2  Assimilation  of  Derived  Density  Fields 
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Figure  2:  Example  of  data  assimilation:  (left)  CHAMP  observations  over  a  30  minute  assimilation  window;  (right) 
assimilated  state  after  combining  CHAMP  observations  with  GITM  forecast. 
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Estimation  of  F10.7  with  CHAMP  obs 


Figure  3:  Left  plot:  daily  averaged  measured  F10.7  (red  line),  ensemble  average  F10.7  from  the  assimilation  of 
CHAMP  neutral  density  observations  (blue  line)  and  ensemble  standard  deviation  for  FI 0.7  (blue  dashed  line).  The 
oscillations  in  the  estimated  FI 0.7  seems  to  follow  the  day-night  change  seen  by  the  CHAMP  satellite.  Right  plot: 
ensemble  average  density  field  for  October  21  2002  at  0800  hours  UTC. 


To  enhance  the  density  estimation  of  GITM  with  data  assimilation,  we  use  derived  density  observations  from  the 
CHAllenging  Minisatellite  Payload  (CHAMP)  [6,7]  and  Gravity  Recovery  and  Climate  Experiment  (GRACE)  [8] 
missions.  Both  mission  satellites  have  GPS  and  accelerometer  data  that  are  processed  to  estimate  a  neutral  density 
field  [9],  Figure  2  shows  an  example  of  assimilating  a  30  minute  observational  window  into  the  GITM  forecast. 
Given  the  scarcity  of  the  data,  it  is  challenging  to  constrain  the  whole  model  field.  To  address  this  issue,  the 
assimilation  is  used  to  estimate  both  model  variables  and  key  model  parameters  that  influence  the  global  evolution 
of  the  model.  To  select  the  most  relevant  model  parameter(s),  a  sensitivity  study  was  performed  a  reveled  that  a  key 
parameter  is  F10.7  index.  This  result  is  consistent  with  the  physical  representation  of  the  F10.7  index  within  the 
model,  since  it  is  being  used  as  a  proxy  for  solar  activity.  The  FI 0.7  index  is  not  strictly  an  independent  input 
parameter  for  the  model,  since  it’s  a  measured  index.  Therefore,  the  assimilation  will  “calibrate”  this  index  by 
adding  a  correction  provided  by  the  data  assimilation.  That  is,  the  assimilation  estimates  d FI0  7  for  the  model  input 
index  in  the  following  way, 

*To.7  =  *10.7  +  <^10.7  (1) 

where  F"'10.7  is  the  model  index  used  by  GITM,  and  F10j  is  the  observed  index.  In  this  formulation,  the  correction  to 
F10.7  is  actually  taking  into  account  model  bias  and  trying  to  correct  it  using  data  assimilation. 


Figure  3  shows  the  daily  averaged  measured  F10.7  value  (red  line),  the  ensemble  average  of  the  estimated  F10.7 
(blue  line)  using  data  assimilation  with  CHAMP  data,  as  well  as  the  ensemble  standard  deviation  of  F10.7.  All  index 
values  are  valid  for  October  21-31  2002.  The  estimated  FI 0.7  value  (blue  line)  is  oscillating  in  accordance  with  the 
day-night  position  of  the  CHAMP  satellite,  indicating  that  the  assimilation  is  correcting  an 
overestimation/underestimation,  suffered  by  the  GITM  model,  through  the  F10.7  index.  The  left  plot  shows  the 
ensemble  average  neutral  density  estimation  for  October  21  2002  at  0800  hours  UTC. 

The  assimilation  results  indicate  that  using  data  assimilation  can  reduce  the  forecast  error  of  the  model.  Furthermore, 
the  model  bias  can  be  corrected  through  the  calibration  of  key  model  parameters  through  data  assimilation.  Other 
key  parameters  that  influence  other  fields  and  properties  of  the  mode,  such  as  temperature,  composition,  etc.,  will  be 
explored  and  included  for  calibration  in  the  data  assimilation.  Initial  condition  and  boundary  conditions  are  also  a 
concern  for  the  model  simulation,  and  will  be  addressed  in  the  assimilation  scheme  as  well. 


4.  DRAG  COEFFICIENT  MODELING 


Drag  is  the  largest  source  of  uncertainty  in  the  orbital  trajectories  of  satellites  in  low  Earth  orbit  (LEO).  The 
acceleration  on  a  satellite  due  to  drag,  an  ,  is  defined  by  [10] 


1  „  A  2 
ao  =  -pCD—v;el 

2  m 


Vrel 


\Vrel 


where  X  is  the  local  atmospheric  mass  density,  Co  is  the  satellite  drag  coefficient,  A  is  the  projected  area  of  the 
satellite  normal  to  the  velocity  vector,  m  is  the  satellite  mass,  and  vrei  is  the  relative  velocity  between  the  satellite  and 
the  co-rotating  atmosphere.  For  artificial  satellites,  the  primary  source  of  drag  acceleration  uncertainty  stem  from 
inadequate  knowledge  of  X  and  Co- 

Atmospheric  mass  densities  are  often  inferred  from  the  orbital  decay  of  satellites  by  using  fixed  or  fitted  drag 
coefficients  [11,12];  however,  the  use  of  such  drag  coefficients  introduces  biases  into  the  inferred  mass  densities 
[1].  These  biases  can  be  partially  removed  by  using  physical  drag  coefficients  computed  based  on  the  momentum 
exchange  between  atmospheric  particles  and  the  satellite  surface.  The  vast  majority  of  LEO  satellites  orbit  in  free 
molecular  flow  where  intermolecular  collisions  between  atmospheric  particles  can  be  neglected.  In  such  conditions, 
the  drag  coefficients  for  satellites  with  simple  convex  geometries  have  closed- form  solutions  [13-15].  A  numerical 
method,  such  as  Test  Particle  Monte  Carlo  (TPMC),  is  required  to  accurately  calculate  drag  coefficients  for  satellites 
with  complicated  concave  geometries.  At  lower  altitudes,  where  the  free  molecular  flow  assumption  breaks  down,  a 
rarefied  gas  dynamic  technique  that  accounts  for  intermolecular  collisions,  such  as  Direct  Simulation  Monte  Carlo 
(DSMC),  is  required  [16]. 

Physical  drag  coefficients  are  dependent  on  a  variety  of  atmospheric  and  satellite  properties,  including: 
atmospheric  translational  temperature,  satellite  surface  temperature,  relative  velocity  between  the  satellite  and  the 
co-rotating  atmosphere,  atmospheric  composition,  satellite  surface  composition,  and,  most  importantly,  the  gas- 
surface  interaction  (GSI)  [13,15].  The  GSI  is  generally  controlled  by  one  or  more  momentum  or  energy 
accommodation  coefficients,  depending  on  the  GSI  model.  One  of  the  simplest  GSI  models,  Maxwell’s  model,  is 
not  dependent  on  accommodation  coefficients.  Instead,  Maxwell’s  model  splits  the  population  of  reflected  particles 
into  a  fraction,  e,  that  is  reflected  specularly,  while  the  remaining  portion,  1-e,  are  diffusely  reflected.  Recently,  we 
have  recently  shown  in  [17]  that  Maxwell’s  model  is  incapable  of  matching  fitted  drag  coefficients  as  a  function  of 
altitude.  Two  more  sophisticated  GSI  models  are  diffuse  reflection  with  incomplete  accommodation  (DRIA)  [18] 
and  the  Cercignani-Lampis-Lord  (CLL)  model  [19].  The  DRIA  model  has  been  applied  in  satellite  drag  coefficient 
modeling  for  nearly  50  years;  however,  the  CLL  model  was  only  recently  applied  to  satellite  drag  coefficient 
modeling  for  first  time  by  our  group  [17]. 

We  have  developed  closed-form  solutions  [17]  for  simple  convex  geometries  such  as  a  sphere  and  flat  plate  for 
use  with  the  CLL  GSI  model.  Closed-form  solutions  were  developed  by  fitting  analytic  expressions  (modified  from 
the  original  Schaaf  and  Chambre  [13]  solutions)  to  DSMC  simulations  with  NASA’s  DAC  [20]  that  were  sampled 
from  the  global  parameter  space  using  Latin  Hypercube  sampling.  These  CLL  closed-form  solutions  fit  the  DAC 
simulations  within  -0.5%  of  the  global  parameter  space. 

Decades  of  work  by  Moe  et  al.  [21-25]  has  shown  that  satellite  surfaces  are  likely  covered  by  varying  levels  of 
atomic  oxygen.  The  atomic  oxygen  adsorbate  changes  the  nature  of  the  GSI,  specifically,  the  effective  energy 
accommodation  coefficient,  a  in  the  DRIA  model,  when  compared  to  particles  interacting  directly  with  a  clean 
spacecraft  surface.  Pilinski  et  al.  [26]  was  the  first  to  fit  the  variation  of  a  as  a  function  of  altitude  with  a  Langmuir 
isotherm  dependent  on  the  partial  pressure  of  atomic  oxygen,  P0\  however,  their  original  model  had  several 
deficiencies,  including  a  going  to  zero  as  Po  goes  zero.  Walker  et  al.  [17]  and  Pilinski  et  al.  [27]  have  improved 
upon  this  deficiency  by  computing  a  as  a  weighted  sum  of  unit  accommodation  ( aajs  -  7)due  to  the  adsorbate 
covered  surface)  and  the  energy  accommodation  coefficient  of  the  clean  satellite  surface,  aa(js,  based  on  Goodman’s 
formula  [28] 

a  =  (l-9)ocsrf  +  0aads 
2.4/7 

(Y  —  - - - 

srf  (X  +  pf 

We  compared  CLL  and  DRIA  GSI  models  [17]  while  including  the  effects  of  atomic  oxygen  adsorption  and 
found  that  both  models  can  match  fitted  drag  coefficients  [29]  equally  well  below  -500  km  altitude.  They  also  found 
larger  (greater  than  -20%)  variations  in  the  drag  coefficients  of  simple  geometries  between  solar  minimum  and  solar 


maximum,  confirming  previous  results  from  Moe  et  al.  [22].  We  investigated  also  in  [17]  the  differences  in  drag 
coefficients  computed  using  NRLMSISE-00  [1]  and  Global  Ionosphere-Thermosphere  Model  (GITM)  [2] 
atmospheric  properties.  At  solar  maximum,  minor  drag  coefficient  discrepancies  of  ~3%  were  found;  however, 
much  larger  differences  (up  to  ~1 1%)  were  found  during  solar  minimum. 

More  recently,  we  developed  a  response  surface  model  (RSM)  technique  [30]  for  drag  coefficient  modeling. 

The  method  was  validated  for  a  sphere  and  then  extended  to  the  more  realistic  case  of  the  GRACE  satellite. 
Comparison  of  the  original  TPMC  training  simulations  and  the  RSM  predictions  show  errors  of  -0.25%  for  the 
sphere  and  -0.7%  for  GRACE.  We  compared  in  [31]  Langmuir,  Temkin,  and  Freundlich  adsorption  models  and 
showed  that  the  Temkin  and  Freundlich  models  match  fitted  drag  coefficients  better  and  both  low  and  high  altitudes. 
Furthermore,  both  adsorption  models  alleviate  some  of  the  physical  deficiencies  present  in  the  Langmuir  model, 
such  as  constant  adsorption  energy  and  monolayer  adsorption. 

5.  ATMOSPHERIC  DENSITY  RECONSTRUCTION 

The  ground-based  tracking  observations  are  used  to  estimate  the  orbital  state  and  drag  ballistic  coefficient  of  a 
number  of  satellites.  By  analyzing  the  change  in  the  satellites’  orbits  over  time,  one  can  estimate  the  atmospheric 
neutral  density,  in  the  form  of  corrections  to  an  assumed  density  model.  This  approach  is  often  called  a  Dynamic 
Calibration  of  the  Atmosphere  (DC A)  in  the  literature  [10].  This  “nowcast”  density  estimate  can  then  be  fed  into  the 
data  assimilation  described  above  to  provide  better  physics-based  density  forecasts. 

We  have  developed  a  new  DCA  method  called  Satellite  Orbit  Tomography,  which  was  originally  inspired  by  X-ray 
computed  tomography.  Here,  rather  than  using  the  decay  in  X-ray  intensity,  we  use  the  decay  in  orbital  specific 
mechanical  energy  (£).  Shoemaker  et  al  [32]  describes  the  mathematical  formulation,  and  Shoemaker  et  al  [33]  gives 
simulation  results  of  an  operationally  realistic  scenario.  The  method  is  outlined  as  follows: 

•  Identify  a  set  of  target  satellites  for  tracking,  nominally  those  that  are  inactive  (e.g.  debris,  rocket  bodies). 
Track  each  target  over  a  span  of  time  (>days)  to  build  up  estimates  of  the  position,  velocity,  and  drag 
ballistic  coefficient  ((3).  In  the  simulations  considered  to  date,  we  have  used  a  Constrained  Admissible 
Region  Multiple  Hypothesis  Filter  (CAR-MHF)  [34-36]  to  estimate  these  states  based  on  measured  angles 
and  angle-rates. 

•  Remove  bias  in  the  estimated  [3  introduced  by  global  errors  in  the  assumed  density  model  used  in  the  CAR- 
MHF.  This  is  done  by  including  at  least  some  tracking  targets  that  have  fairly  well-modeled  [3,  and 
comparing  their  modeled  [3  with  that  estimated  by  the  CAR-MHF.  This  approach  also  leverages  the  drag 
coefficient  modeling  described  above:  given  the  published  information  on  a  satellite’s  shape  and  mass,  we 
are  able  to  sufficiently  model  [3. 

•  Using  a  given  satellite’s  estimated  position  and  velocity  at  one  time  (//),  and  then  again  at  a  later  time  ( t2 ), 
calculate  the  decay  in  £,  using  the  osculating  orbit  states  at  those  times.  The  time  span  At  =  t2-  tj  should  be 
long  enough  to  observe  the  decay  signal  above  the  measurement  noise,  yet  short  enough  to  recover  some 
time-resolved  information  in  the  density  model.  Our  simulations  have  used  At  of  48  hours,  but  a  lower  limit 
of  24  hours  is  feasible  given  the  expected  accuracy  of  our  ground-based  tracking  system  and  the  likely 
revisit  rate  (-1  to  2  passes  per  night)  for  our  mostly  LEO  targets. 

•  Using  the  measured  decay  in  £,  from  the  set  of  targets,  solve  for  a  spatially-resolved  scalar  correction  ( s )  to 
the  assumed  density  model.  The  correction  factor  .?  is  defined  in  a  grid;  we  have  used  grids  spanning  300  to 
500  km  altitude,  with  100  km  altitude  spacing  and  20  deg  spacing  in  latitude  and  longitude.  In  general,  the 
problem  is  underdetermined  (there  are  more  grid  elements  than  target  satellites)  and  ill-posed  (most 
satellites  do  not  pass  through  each  grid  element).  Thus,  we  use  Tikhonov  regularization  to  stabilize  the 
solution,  with  a  spatial  smoothness  constraint  on  the  solved-for  s  field. 

Our  most  recent  simulations  use  actual  resident  space  objects  from  the  publically  available  catalog,  and  assume  a 
single  ground  site  at  Fenton  Hill,  NM.  Using  only  40  targets  and  having  suitable  visibilities  from  this  single  ground 
site,  we  are  able  to  reconstruct  the  time-averaged,  yet  spatially  resolved,  density  field  over  48  hours  to  within 
approximately  10%.  These  results  also  assume  reasonable  orbit  estimation  errors,  and  that  each  satellite’s  assumed 
ballistic  coefficient  [3  will  have  zero-mean  error  with  1-a  standard  deviation  of  10%.  The  satellite  orbit  tomography 
approach  has  some  practical  advantages  over  typical  weighted  least-squares  approaches,  such  as  allowing  easy 
density  model  specification  i.e.  not  requiring  Jacobian  matrices  that  describe  the  sensitivity  of  the  density  model 
dynamics  to  the  states. 


6.  MACHINE  LEARNING  FOR  FAST  DENSITY  FORECASTING 


Given  the  significant  computational  expense  of  a  full  physics-based  simulation  of  the  spatial  variation  of  the 
atmospheric  density  profile,  there  is  a  need  for  a  rapidly  computable  model  approximating  this  profile.  Our  goal  is  to 
apply  machine-learning  methods  to  construct  such  a  model,  enabling  estimation  of  the  density  profile  based  on 
current  measurements  of  physical  parameters  such  as  solar  activity  and  terrestrial  magnetic  field.  Our  initial  attempts 
to  learn  a  model  based  on  historical  point-measurements  of  atmospheric  density  from  instrumented  satellites  (e.g. 
CHAMP  and  GRACE)  showed  that  the  spatial  sampling  density  of  available  data  is  entirely  inadequate  for  such  a 
task.  Our  current  approach,  therefore,  is  to  construct  a  computationally  cheap  surrogate  model  learned  from 
simulations  produced  by  the  GITM  code;  this  work  is  still  in  progress.  Our  initial  task  has  been  the  construction  of  a 
low-dimensional  parameterization  of  the  density  profile;  given  the  difficulty  of  directly  estimating  a  full,  high¬ 
dimensional  representation  of  the  profile,  a  regression  model  will  be  constructed  to  estimate  the  low-dimensional 
parameterization  from  the  available  physical  measurements.  While  spherical  harmonics  are  a  very  widely  used  basis 
on  which  to  represent  these  profiles,  we  are  currently  investigating  the  use  of  a  basis  derived  from  our  training  data 
(i.e.  a  set  of  GITM  simulations)  via  Principal  Component  Analysis  (PCA). 

7.  UNCERTAINTY  QUANTIFICATION 


Our  goal  is  to  compute  collision  probabilities  that  account  for  uncertainty  at  each  stage  of  our  modeling  procedure. 
To  be  accurate,  a  collision  probability  must  necessarily  incorporate  uncertainty  from  observations,  density  modeling 
and  forecasting,  drag  estimation,  and  other  sources.  In  order  to  combine  uncertainty  across  these  sources,  we  favor 
Monte  Carlo-based  methods. 


Figure  4:  (left)  Gaussian  approximations  to  Monte  Carlo  samples  for  a  simulated  example,  (middle)  Mixture  of 
Gaussian  approximation  to  Monte  Carlo  sample  from  a  simulated  example,  (right)  The  yellow  and  green  points 
show  100  importance  sampled  points  using  the  two-stage  procedure.  This  small  sample  produced  a  very  accurate 
estimate  of  the  probability  when  compared  to  a  brute  force  simulation  of  9,000,000  pairs. 

The  recent  FORTE/Meteor  close  approach  on  22  June  2013  provided  an  opportunity  to  test  a  basic  Monte  Carlo 
strategy.  First,  we  sampled  ten  GITM  densities  from  the  distribution  of  possible  densities.  Next,  we  sampled  100 
satellite  state  vectors  from  the  distribution  induced  by  observations  of  the  two  satellites  made  before  the  close 
approach.  These  state  vector  samples  were  randomly  paired  with  a  sampled  GITM  density  and  the  pair  was 
propagated  to  close  approach.  At  the  time  of  close  approach,  we  approximated  the  distribution  of  the  positions  with 
a  simple  Gaussian  and  used  this  to  calculate  the  probability  of  collision.  This  method  required  relatively  few 
simulations  and  incorporated  two  important  sources  of  uncertainty.  More  uncertainty  sources  could  be  incorporated 
very  easily  (we  calculated  a  best  fit  drag  coefficient,  but  weren’t  able  to  compute  uncertainty  at  the  time),  although 
this  would  likely  require  a  larger  set  of  simulations.  Figure  4  (left)  shows  a  fit  based  on  a  similar  simulated  close 
approach. 

A  simple  extension  of  the  previous  method  is  the  use  of  a  mixture  of  Gaussian  model  for  the  positions  at  close 
approach.  This  collision  probability  would  be  more  accurate  since  the  mixture  model  would  be  a  better 
representation  of  the  final  position  uncertainty.  The  calculation  is  a  straightforward  extension  of  the  single  Gaussian 
case.  One  drawback  is  that  this  method  would  require  a  larger  set  of  simulations  in  order  to  accurately  estimate  the 
mixture  model.  Figure  4  (middle)  shows  the  mixture  of  Gaussian  fit  to  the  same  simulation  as  the  figure  on  the  left. 
Notice  that  this  fit  is  better  at  capturing  the  tails  and  the  overly  dense  central  region. 


We  are  also  working  on  an  importance  sampling  methodology  for  computing  collision  probabilities  with  reduced 
simulation  costs.  Importance  sampling  is  a  Monte  Carlo  method  that  draws  a  biased  sample  that  is  weighted  to 
produce  the  correct  expectations.  Particularly  in  the  case  of  small  probabilities,  this  technique  can  be  used  to  reduce 
the  variance  of  an  estimate,  and  therefore  the  number  of  samples  required.  We  developed  a  two-stage  procedure.  In 
the  first  stage  a  raw  Monte  Carlo  sample  is  drawn  from  all  of  the  uncertainty  sources.  The  results  of  this  stage  are 
analyzed  to  find  regions  of  the  original  distribution  that  produce  close  approaches  and  the  second  stage  samples 
directly  from  these  regions.  Figure  4  (right)  shows  the  results  of  the  methodology  applied  to  the  simulations  shown 
in  figures  on  the  left  and  middle.  The  second  stage  used  100  simulations  to  accurately  estimate  the  collision 
probability  when  compared  with  a  brute  force  simulation. 

8.  MODEL  INTEGRATION  AND  VISUALIZATIONS 

In  order  to  coherently  integrate  these  diverse  models  (Figure  5),  plus  alternative  models  to  facilitate 
comparisons  of  model  techniques  and  implementations,  we  have  developed  an  integration  tool.  For  this  we  chose  the 
Python  language  for  its  rapid  prototyping  and  its  multilingual  extensibility.  Our  tool,  as  open-source,  pure-Python 
package  called  sysdevel  [37],  handles  three  aspects  of  this  multi-model  synthesis: 

1.  unifying  potentially  disparate  model  data  in  a  scalable  and  malleable  way, 

2.  providing  a  service  framework  for  plotting  and  visualizing  data  in  a  web  browser, 

3.  extending  the  Python  build  system  to  include  external  library  dependencies,  unit  testing,  and  automatic 
documentation  to  ease  distribution  of  collaborators  who  may  use  differing  compute  environments. 

The  sysdevel  build  system  extends  the  built-in  Python  distutils  package  to  recursively  build  sub-packages  that  each 
build  one  of  our  models  with  a  normal  'python  setup .  py  build'  call.  This  then,  for  example,  descends  into 
our  GITM  sub-package  like  a  recursive  'make'  call.  It  locates  GITM's  MPI  and  HDF5  library  dependencies  and  the 
proper  Fortran  compiler,  fetching  and  installing  them  if  any  are  missing  by  utilizing  CMake-style  configuration  files 
in  sysdevel.  Finally  it  creates  a  native  executable  for  use  in  a  cluster.  Those  familiar  with  Python  will  recognize  that 
this  is  well  beyond  the  normal  distutils  build  process. 
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Figure  5:  IMPACT  Flowchart  illustrating  the  coupling  between  different  models  for  atmospheric  density,  drag, 
propagation,  observations,  and  uncertainty  quantification. 


Within  the  IMPACT  source  tree,  one  of  our  sub-packages  is  labeled  'website'  and  this  handles  the  user  interface  over 
the  web.  Sysdevel  provides  a  Javascript  and  PHP  framework  specifically  for  simulation  configuration  and  results 
viewing.  Using  a  plumbing  metaphor  (through  jsPlumb  [38]),  the  user  graphically  configures  the  simulation 
processing  pipeline  and  chooses  data  plots  of  interest.  These  interactive  results  plots  are  displayed  as  soon  as  data  is 
available.  To  minimize  computation,  we  store  intermediate  results,  so  data  from  a  duplicate  configuration  is 
available  immediately  unless  caching  is  overridden.  Because  the  framework  communicates  with  the  simulation 
using  JSON  over  WebSockets,  we  were  able  to  extend  our  UI  (outside  of  sysdevel)  to  include  a  3D  visualization 
(built  with  Three.js  [39])  of  RSO  tracks  around  the  Earth,  as  seen  in  Figure  6,  with  relative  ease. 

By  far  the  most  interesting  portion  of  sysdevel  is  its  model  integration  facility.  With  this,  we  build  a  script  that  acts 
as  a  server  to  the  web  interface  or  a  direct  command-line  interface.  Sysdevel  combines  our  simulation  models 
together  using  the  Model-View-Controller  pattern,  which  consists  of  data  models  (in  our  case,  RSO  objects),  data 
controllers  that  manipulate  those  objects  (our  modeling  processes),  and  data  views  that  represent  direct  hooks  to  our 
plotting  and  visualization.  Note  the  conflicting  semantics  for  the  word  'model'  above  -  in  the  context  of  this  well- 
known  software  pattern;  we  will  substitute  the  word  'object'.  These  objects  are  the  cores  of  our  model  synthesis 
strategy.  As  our  development  process  of  IMPACT  iterates  from  1 -to- 1 -conjunction  analysis,  expanding  to  all-to-all 
RSO  collision  detection  is  challenging  but  the  additional  complexities  appear  feasible.  Our  data  object  abstraction 
supports  multiple  data  storage  backends  to  allow  for  scalability  of  the  analysis.  Currently  this  consists  of  structured 
Hierarchical  Data  Format  (HDF)  fdes,  but  could  also  utilize  a  relational  database.  As  we  expand  to  cover  more  and 
more  RSOs  and  over  wider  time  scales,  our  storage  performance  needs  grow  drastically.  We  are  currently 
considering  graph  databases  to  meet  that  need,  and  our  abstraction  layer  makes  that  adjustment  possible. 

We  are  not  only  concerned  with  scalability,  but  also  malleability.  We  already  have  integrated  some  alternate  models 
(such  as  MSIS  for  atmospheric  density),  but  we  want  to  easily  include  other  alternatives  throughout  the  pipeline 
without  extensive  integration  work.  To  that  end,  sysdevel  data  objects  are  self-describing  using  built-in  Python 
idioms.  This  feature  allows  us  to  simply  alter  our  data  object  definition  (by  creating  a  new  sub-class  of  the  original 
Python  data  object)  to  also  conform  to  the  domain  ontology  of  the  new  model,  and  the  user  is  ready  to  go.  Inside 
sysdevel,  there  is  a  great  deal  more  complexity  to  map  that  change  to  the  storage  backend  (hence  our  initial 
preference  for  HDF,  which  simplifies  this  mapping). 

Through  these  features  of  scalable  and  malleable  data  unification,  simulation  pipeline  configuration  and  data 
visualization,  and  a  comprehensive  build  system,  our  sysdevel  integration  package  not  only  serves  the  needs  of  our 
IMPACT  project,  but  also  provides  a  general  tool  for  other  multi-model  simulations  that  would  otherwise  require 
extensive  effort  to  tie  together. 


Figure  6:  Screen  snapshot  of  the  visualization  interface  of  the  IMPACT  framework.  Satellite  position  is  represented 

with  a  particle  cloud  indicating  the  positional  uncertainty. 


9.  SUMMARY 


The  IMPACT  project,  funded  by  Los  Alamos  National  Laboratory  with  internal  research  and  development  funding, 
has  developed  an  open  source  and  flexible  research  framework  for  analyzing  satellite  conjunctions  with  modem 
physics-based  atmospheric  density  models,  drag  coefficient  analysis,  a  new  atmospheric  density  reconstruction 
method,  machine  learning,  and  system  based  uncertainty  quantification.  We  presented  an  overview  of  this  project 
including  some  of  the  novel  approaches  applied  to  satellite  conjunction  analysis.  We  welcome  future  collaborations 
using  the  IMPACT  framework. 
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